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ABSTRACT 

We report on the first study of the angular distribution of energetic par- 
ticles and radiation generated in relativistic collisionless electron-positron pair 
plasma reconnection, using two-dimensional particle-in-cell simulations. We dis- 
cover a strong anisotropy of the particles accelerated by reconnection and the 
associated strong beaming of their radiation. The focusing of particles and ra- 
diation increases with their energy; in this sense, this "kinetic beaming" effect 
differs fundamentally from the relativistic Doppler beaming usually invoked in 
high-energy astrophysics, in which all photons are focused and boosted achro- 
matically. We also present, for the first time, the modeling of the synchrotron 
emission as seen by an external observer during the reconnection process. The 
expected lightcurves comprise several bright symmetric sub-flares emitted by the 
energetic beam of particles sweeping across the line of sight intermittently, and 
exhibit super-fast time variability as short as about one tenth of the system light- 
crossing time. The concentration of the energetic particles into compact regions 
inside magnetic islands and particle anisotropy explain the rapid variability. This 
radiative signature of reconnection can account for the brightness and variability 
of the gamma-ray flares in the Crab Nebula and in blazars. 

Subject headings: Acceleration of particles — Magnetic reconnection — Radia- 
tion mechanisms: non-thermal — ISM: individual (Crab Nebula) — Galaxies: 
active 
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Introduction 



The high-energy radiation from numerous astrophysical objects, including active galac- 
tic nuclei, pulsar wind nebulae, and gamma-ray bursts, is emitted by particles accelerated 
to relativistic speeds. Magnetic reconnection is one of the main mechanisms thought to 
accelerate par ticles, by converting mag netic energy into particle kinetic energy (see, e.g., 
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tailed picture of the particle energy spectrum. However, by itself, the energy spectrum lacks 
information critical to the determination of radiation generated by the plasma — namely, 
the angular distribution of the velocities of accelerated particles. Because ultra-relativistic 
particles radiate in a narrow cone along their direction of motion, any anisotropy of the 
energetic particles translates directly into the anisotropy of their emission (e.g., synchrotron 
or inverse Compton). The beaming of the radiation drastically affects how we infer, from 
observations, the physical conditions of the emitting region (e.g., size, overall energetics and 
dynamics) and statistical properties of flaring astrophysical objects. 

In this Letter, we report on the first detailed analysis of the angular distribution (in addi- 
tion to the energy and spatial distributions) of particles accelerated in collisionless relativistic 
pair reconnection, using PIC simulations. In Section [2] we describe the simulation setup. In 
Section [31 we present our results and report on the discovery of a "kinetic beaming", i.e., 
a strong energy- dependent anisotropy of the particles and their radiation. We also present, 
for the first time, the modeling of the high-energy radiation spectrum and lightcurve as seen 
by a distant observer, and predict extremely rapid time-variability (much shorter than the 
light-crossing time of the system). In Section HJ we briefly discuss the general implications 
of our findings in the astrophysical context of flaring high-energy gamma-ray sources like the 
Crab Nebula and blazars. 



2. PIC simulation setup 



We performed two-dimensional numerical simulations of collisionless relat ivistic pair 
plasm a reconnection using the explicit electromagnetic PIC capabilities of vorpal ( iNieter &; Cary 
2004T). The initial setup adopted here is standard in reconnection simulations (see, e.g., 



Zenitani fc Hoshino! 1200 ll ). It consists of a rectangular box of size L Y x L y with two anti- 



parallel relativistic Harris current layers (IKirk &; Skjaeraasenl 120031 ) and double periodic 
boundary conditions. In the following, we will focus on the dynamics of the bottom layer 
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only, i.e., the bottom half of the simulation domain. The reconnecting magnetic field is 
B x = Bq tax\h(y / S) , where -Bo is the upstream field and 5 is the initial layer thickness. There 
is no guide field, B z = 0. The simulation has a resolution of 4.6 grid cells per 5 ~ 0.8p c , 
where p c = m e c 2 /eB is the non-relativistic electron Larmor radius, m e is the electron rest 
mass, e the elementary electric charge, and c the speed of light. 

The initial particle density (electrons and positrons together in the laboratory frame) is 
n = ridrift+^O; where n^ft = cosh~ 2 (?//5) is the density of electrons and positrons drifting 
in opposite directions (in the indirection) at a velocity ^drift/ c = 0.6 and located in the 
layer, and uq = 0.042rido is a uniform and isotropic background density. Both populations are 
distributed according to a relativistic Maxwellian of temperature kT' drih = 0.3m e c 2 defined in 
the drifting particles co-moving frame, and kT^ g = 0.15m e c 2 in the laboratory frame, where 
k is the Boltzmann constant. The electron skin depth in the layer is d c = c/u pc ~ 0.7 p c , 
where oo pe = y 'Airridoe 2 / 'm c is the plasma frequency defined with the drifting electrons and 
positrons. The upstream plasma beta parameter is (3 = 8irnokTb g / B 2 ps 2.6 x 10~ 2 (the 
magnetization sigma parameter is a = B 2 1 ^imom P c 2 ~ 11.5). Radiative energy losses are 



neglected (see iJaroschek fc Hoshinoll2009l for a PIC simulation including radiative drag). 



The unit of length in the x- and ^-directions is p c , and the unit of time is the inverse of 
the nominal cyclotron frequency u~ l = p c /c. In order to initiate the reconnection process, 
we break the initial Harris equilibrium with a tearing-like perturbation in the magnetic flux 
function (Figure dj top panel). We performed a series of simulations with different box 
sizes {L x /p c ,L y /p c ) = (90,90), (180,180), (360,360), and (720,720), and with 16, 64, and 
256 particles per cell. We verified numerical convergence with respect to both the spatial 
resolution (i.e., the number of grid cells per p c ) and the number of particles per cell. The 
results shown below are for a L x = L y = 360p c and 2.7 x 10 8 particles (64 per cell). 



3. Results 



3.1. Overall reconnection dynamics 



Figure [T] shows the three main stages of the reconnection dynamics for the bottom 
layer only. The initial current layer quickly becomes unstable to secondary tearing mode s 



[e.g., IJaroschek et al.ll2004bl ; lLoureiro et al.1 120071 ; iKomissarov et al.l 120071 ; iLiu et al.l 120111 1. 
breaking into a dynamical chain of magnetic islands made of closed magnetic flux loops, 
visible in the intermediate stage at t = Sldu^ 1 ~ 0.9L x /c (middle panel in Figure [TJ), 
when about 44% of the initial magnetic energy has been dissipated. The number of islands 
subsequently decreases as small islands merge into bigger ones. We estimate the time- 
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Fig. 1. — Plasma density map and magnetic field lines (white solid lines) at the initial stage 
tu c = (top panel), at the intermediate stage tu c = 319 (middle panel), and close to the 
final stage of the reconnection event tu c = 637 (bottom panel). Only the bottom half of the 
simulation box is shown here. Distances in the x- and y-directions are normalized to the 
non-relativistic electron Larmor radius p c . The density is normalized to the initial drifting 
particle density n d0 . 



-5 - 



averaged dimensionless reconnection rate to be /3 rec = E z /Bq w 0.16 (where E z is the electric 
field in the z-direction) , up to t = QOOu' 1 when most of the field has reconnected. At the 
end of the reconnection event, the system settles in a saturated configuration with a single 
big island, or O-point, and a single X-point (see bottom panel in Figure [TJ t = GSTcj" 1 ) in 
each half of the box. The final state in the upper-half of the simulation domain (not shown 
here) is identical to the bottom-half but symmetric with respect to the center of the box. 
All the magnetic flux ends up around the two main islands, and the separatrices emanating 
from the two X-points connect to each other. During this process, about 55% of the initial 
magnetic energy is converted into particle kinetic energy. The total energy in the system is 
conserved with less than 0.1% of error at the end of the simulation t = 1270a; ( r 1 . 



3.2. Particle energy distribution, particle anisotropy and spatial 

inhomogeneity 

Shortly after the onset of reconnection, a new population of relativistic particles emerges 
from the initial cool thermal distribution. Particles are accelerated along the z-direction at 
X-points by the strong reconnection electric fi eld E 7 .. They are then deflected along the ±x- 



direction by the reconnected magnetic field B y (IZenitani fc Hoshinoll200ll ; ISironi fc Spitkovsky 



201ll ). Figure [2] (upper panel) presents the energy distribution / -f 2 dN/d / ~f of all the positrons 
in the simulation box (blue solid line), where 7 is the particle Lorentz factor, at the interme- 
diate stage of reconnection t = 319a;" 1 . The high-energy bump peaking at 7 ~ 8 can be well 
fitted by dN/dj oc 7~ 1//2 exp(— 7/5). We describe this component as a quasi-thermal dis- 
tribution resulting from plasma heating by magnetic dissipation, rather than a non-thermal 
power-law tail. It contains about 49% of the particles and 78% of the kinetic energy. 

The main result of this paper concerns the energetic particles' anisotropy, which we 
examine by calculating the total solid angle within which half of the particles of a given 
energy are contained, fi e ,5o%(7); as a function of the particle energy (FigureEl bottom panel). 
The low-energy particles (7 < 10) remain approximately isotropic, whereas the high-energy 
particles (7 > 20) are focused in a tight beam whose solid angle decreases rapidly with energy, 
roughly as fi e ,5o% oc 7~ 3,5 . The beam's angular size can become as small as f2 ei 5o%/47r < 1% 
of the whole sphere for 7 > 40. 

In addition to the anisotropy, we study the spatial distribution of the energetic particles 
in the computational domain. We quantify the degree of inhomogeneity of the particles by 
computing the fraction of the system surface (L x x L y ) covered by half of the particles, as a 
function their energy, / e ,5o%(7) (Figure El bottom panel). We find a strong energy- dependent 
inhomogeneity: the particles with 7 < 2 fill a large fraction of the box (/ e ,5o% ~ 0.4), mostly 
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Fig. 2. — Top panel: Energy distribution function r y 2 dN / of all the positrons in the 
simulation box (blue solid line) as function of 7 at tu c = 319. The dotted lines are analytical 
fits to the low- (relativistic Maxwellian of temperature kT = 0.15m e c 2 ) and high-energy 
(dN/dj oc 7 -1 / 2 exp(— 7/5)) parts of the energy distribution function. Bottom panel: Solid 
angle normalized by 4n, f2 e;5 o%/47r (red dot-dashed line), and spatial filling factor, / e ,5o% 
(blue solid line), containing half of the positrons in a given energy bin, as functions of 7 at 
tu c = 319. The pink dashed line is a power-law fit to n e ,50%/47r for 7 > 20. The gray bands 
are three particle energy bins (a), (b), and (c) for which the angular and spatial distributions 
of the particles are shown in Figure |3j 
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Fig. 3. — Angular distribution maps (using the Aitoff projection, panels (al)-(cl)), and 
spatial distribution maps (panels (a2)-(c2)) of the positrons in the three energy bands (a)- 
(c) defined in Figure El In the angular maps, the x-axis gives the value for the longitude A 
and the y-axis the latitude </> (see the text for their definitions). In panel (al), the directions 
±x, ±y and ±z are indicated. The subdomains —15° < 4> < +15°, —105° < A < —75° 
labeled "(1)" and +30° < <\> < +60°, +45° < A < +75° labeled "(2)" in panel (cl) shown 
by the white arrows and delimited by white squares are used in Figure H] to compute the 
synchrotron radiation spectrum emitted in these specific directions. Bright/dark colors show 
high/low densities of particles per unit of solid angle (left panels) and per unit of area (right 
panels). 
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outside magnetic islands, while the higher energy particles (2 < 7 < 20) are co ncentrated 



i nto s mall bunches (f e ,5o% < 0.1) inside islands. These results are consistent with iLiu et al. 
(120 111 ) simulations. For 7 > 20, / c .5o% drops abruptly. 

Figure [3] provides a good illustration of the anisotropy and inhomogeneity of the particle 
distribution in the three different energy bins shown in Figure [21 The angular distribution 
is calculated using spherical coordinates in which a radial unit vector has the coordinates 
x = cos sin A, y = sin0, z = cos </> cos A, where A is the longitude and <fi the latitude. After 
the end of reconnection, the energy dependence of anisotropy and inhomogeneity decreases 
because all particles end up circling inside the O-point (Figure (H bottom panel). 



3.3. Synchrotron beam 

Next, we examine the main radiative signatures of the reconnection process, namely, the 
emission spectrum and temporal variability. We first look at the effect of particle anisotropy 
on the observable synchrotron radiation spectrum emitted by the layer. Following the same 
procedure as for the particle energy and angular distributions, we characterize the photon 
distribution emitted by the particles. Our analysis makes four approximations: (1) the 
particles emit pure synchrotron radiation, (2) the plasma does not absorb the radiation 
(optically thin), (3) all the emission is beamed in the direction of motion of the radiating 
particle (valid for 7>1), and (4) synchrotron energy losses and the radiation reaction force 
on the particles are ignored. All the results presented below regarding the calculation of 
radiation are performed after the simulation is completed, in accordance with assumptions 
(2) and (4). 



Using the classical synchrotron spectrum formula (e.g.. iBlumenthal fc Gouldlll97Q[ ). we 
calculate the resulting instantaneous photon spectral energy distribution (SED, i.e., radiative 
power per unit of area) emitted by all the positrons in the box at tu c = 319 (Figure H]). 
Frequencies are normalized to the nominal critical synchrotron frequency v c = 3u c /4tt. The 
overall shape of the SED averaged over all directions (vF v )\s (blue solid line) resembles the 
shape of the particle energy distribution in Figure |2j The spectral peak coincides with the 
typical synchrotron photon frequency of the bulk of energetic particles (7 ~ 10), i.e. vjv c ~ 
7 2 = 100. Below the peak [yjv c < 100), the spectrum can be well fitted by a single power law 
of index ~ +0.6. The cool initial distribution of particles (with kT = 0.15m e c 2 ) is responsible 
for the slight flux excess at low frequencies {yjv c < 10). The most energetic particles (7 > 10) 
radiate above vjv c = 100 and form a soft power-law- like component of index ~ — 0.7 between 



200 < vjv c < 2000, followed by a sharp cut-off (see also I Jaroschek et al.ll2004al for a similar 
calculation). 
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Fig. 4. — Instantaneous spectral energy distribution emitted by all the positrons in the 
bottom-half of the simulation box via synchrotron radiation vF v averaged over all directions 
as a function of the reduced frequency v jv c (blue solid line), where v c = 3u; c /47r, at tu c = 319. 
For comparison, the green dashed lines show the spectral energy distributions emitted by 
the particles contained in the solid angle domains (1) and (2) defined in Figure [31 panel (d). 
This figure shows also the solid angle containing half of the photons in a given frequency 
bin, f2 p h,5o%, normalized by 4ir, as a function of vjv c (red dot-dashed line). 
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The anisotropy of the high-energy particles translates directly into the anisotropy of 
radiation. We compute the angular distribution of the emission using the same measures as 
for the particles (see Section l3T2~j) . namely the solid angle within which half of the photons 
are contained, fi p h,5o%; as a function of v jv c (Figure HI red dot-dashed line). As expected, we 
find that the emitted flux displays a strong frequency-dependent anisotropy, very much like 
the particles, although the transition from the isotropic to the highly anisotropic regime, 
roughly at vjv c = 100, is more gradual for photons. The solid angle of the radiation 
beam decreases with frequency approximately as n phj50 % oc (i//V c )~°- 75 . The high-energy 
photons [yjv c > 100) are concentrated in a small solid angle fiph,5o%/47r < 0.10. The 
angular distribution maps (similar to Figure [3J, not shown here) indicate that the high- 
energy radiation is strongly beamed towards the ±x-directions at tu c = 319, although the 
beam is changing direction restlessly within the plane of the layer during reconnection. 

To illustrate the significance of beaming, we present the spectrum of photons (uF„)m 
emitted in the direction of the most energetic particles, e.g., around the — x-direction (—15° < 
<p < +15°, —105° < A < —75° corresponding to Afim = 0.27 sr, see Figure[3j domain "(1)" in 
panel (d)), and compare it with {uF u )^ (FigureS]). The spectrum (uF v )^ is notably harder 
than {vF v )- lso at all frequencies. The beaming of the most energetic particles concentrates 
their synchrotron radiation into a small solid angle, yielding a flux more than an order of 
magnitude greater than the isotropic flux at the same frequency. In contrast, the observed 
high-energy emission is strongly suppressed in other directions as, for instance, in the solid 
angle domain (2) (Figures EH])- The results are qualitatively identical for particles radiating 
predominantly via inverse Compton scattering, because target photons are scattered and 
focused in the direction of motion of the particles. 



3.4. High-energy radiation lightcurve 

Consider an observer at infinity looking in the same direction during the entire reconnec- 
tion event. What would be the high-energy radiation flux seen by the observer as a function 
of time? To calculate the lightcurve, we compute the flux received by the observer taking 
into account the time delay due to the finite time of propagation of the radiation through 
the box. As an example, we consider an observer looking in the directions ±x (A = ±90°, 
= 0°). We sum the contributions from all the particles going in the direction delimited by 



1 Although this is not the case here, if the particle beam solid angle ri .so% were smaller than ~ I/7 2 , 
then the angular size of the radiation beam would be controlled by the opening angle of the synchrotron 
beam of a single particle, i.e., f2 p h,5o% ~ 1/7 2 - 
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Fig. 5. — High-energy synchrotron flux above vjv c = 100 as a function of time (in units of 
the light-crossing time of the simulation box L x / c) seen by an observer located at infinity, 
looking in the — x (red solid line) and +£-directions (blue dotted line) within AQ = 0.03 sr. 
The total lightcurve averaged over all directions is shown for comparison (green dashed line 
labeled "iso."). 
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the finite but small solid angle domain AQ± X = 0.03 sr centered around the ±x-directions. 
Figure [5] gives the observed photon flux integrated above vjv c = 100 as a function of time. 

We find that reconnection generates bright sub-flares on timescales of order one tenth 
the light-crossing time of the system (L x /c). The amplitude of the spikes increases with 
the observed radiation frequency. The short time-variability is due to the bunching of the 
high-energy particles into small volumes inside the magnetic islands moving away from the 
X-points along the ±x-directions, and particle anisotropy. The high-energy beam of particles 
sweeps across the line of sight intermittently and generates each bright spike of the lightcurve 
with nearly symmetric profile (i.e., the rising time is of order the decaying time). The 
intense sub-flares are smoothed out if one considers the total flux averaged over all directions 
(Figure El green dashed line), demonstrating that they are caused by a geometric effect 
(sweeping beam) rather than an intrinsic change in the acceleration mechanism. At the end 
of the reconnection process, even the high-energy variability decays due to the isotropization 
of particles at the O-point. 



4. Astrophysical implications 



The anisotropy of the particle distribution function discovered in this study leads to a 
strong beaming of the radiation emitted during a reconnection event. This "kinetic beam- 
ing" is energy-dependent, i.e., the collimation of particles and radiation increases with their 
energy. Kinetic beaming differs from the relativistic Doppler beaming usually invoked in 
high-energy astrophysics (IReesI Il966h : Doppler beaming is caused by the bulk motion of 
a plasma emitting isotropically in its rest frame; in contrast to kinetic beaming, Doppler 
beaming focuses and boosts all photons by the same factor regardless of their energies. 
This fundamental difference provides a way to discriminate observationally between these 
two beaming mechanisms. In addition, we expect rapid variability of the observed flux 
much shorter than the light-crossing time of the system with nearly symmetric burst profiles 
(particle bunching and sweeping beam). This situation is often encountered in high-energy 
astrophysics, in objects such as, e.g., active galactic nucleus jets, or gamma-ray bursts. 



The discovery of gamma-ray flares in the Crab Nebula (ITavani et al.ll2011t lAbdo et al. 



20111) is a good example, because the shortest detected variability timescale of a few hours 
( iBalbo et al.ll201ll ; iBuehler et al.ll2012l ) may be much shorter than the light-crossing time of 
the flaring region (days to weeks). The nearly symmetric shape of the observed sub-flares 
suggests that the rapid variability is due to a geometric effect. This is consistent with our 
findi ngs, supporting the magnetic reconnection scenario for the origin of the flares in the neb- 
ula (jUzdensky et al .1120 lit ICerutti et al.ll2012f ). in which PeV particles are accelerated and 
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focused in a thin fan beam in the layer. In addition, the super-fast variability (variability 



timescales much shorter than the light-crossing time of the super massive black ho 



served at TeV energy gamma rays in a few blazars (PKS 2155 — 304 (lAharonian et al. 



e) ob 



2003) 



Mrk 501 jAlbert et al.l l2007h or more recently in PKS 1222 + 216 jAleksic et alJboilh ) 



is 



difficult to explain unless one invokes extr eme jet bulk Lorentz factors T > 50 (see, e.g., 
Henri &: Sauee l200d : iBegelman et~alll2008l ). High-energy particle anisotropy and inhomo- 
geneity generated by magnetic reconnection in the comoving frame can alleviate the severe 
constrai nts on the energetics an d collimation of the relativistic jet inferred from TeV obser- 
vations (INalewajko et al.ll2012l ). Finally, the beaming of the high-energy radiation is also 
important for the interpretation of flare statistics. Gamma-ray flares could occur repeatedly 
but we detect only those with emission pointing toward us. 
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